module credmod
	use basicmodule
	implicit none
	
	type tsttype
		real	:: cv_zt(2)
		real, allocatable	:: lam(:), lam_stage1(:), thlist0(:,:), thlist0_stage1(:,:)
		integer	:: n0active, n0active_stage1,stage
	end type
	
	type(tsttype)	:: tst1, tst5, tst10

	real		:: cv_zt(2)

	contains
	
	subroutine settst(tst)
		type(tsttype)	:: tst
		tst%cv_zt=cv_zt
		tst%lam=lam
		tst%lam_stage1=lam_stage1
		tst%thlist0=thlist0
		tst%thlist0_stage1=thlist0_stage1
		tst%n0active=n0active
		tst%n0active_stage1=n0active_stage1
		tst%stage=stage
	end subroutine
	
	subroutine setcvs
		cv_zt(1)=gausscdfinv_s(1-level/2)
		cv_zt(2)=tin(1-level/2,80+10*log(level))
		print *,"cvs", cv_zt
	end subroutine
	
	function switchind(y) result(val)
		real	:: y(1:k), val
		val=max(0.0,y(1)-swsmall)
		if(y(k)<=0) return
		val=min(val,max(0.0,y(1)/y(k)-swflat))
	end function
	
	function rej_stage1(y1,y2,f1,tst) result(val)
		type(tsttype)	:: tst
		real	:: y1(0:k),y2(0:k),f1
		logical :: val
		real	:: h,ma,var,Zsum,msx(3)
		integer	:: j

		Zsum=y1(0)-sum(y2(0:k))
		var=1+sum(y2(1:k)**2)
		h=0
		do j=1,tst%n0active_stage1
			msx=tst%thlist0_stage1(:,j)
			ma=getmeanadj(msx,y1(k))
			h=h+tst%lam_stage1(j)*getdens_base(msx,y1(1:k))*exp(-0.5*(Zsum+ma)**2/var)/sqrt(var)
		enddo
		val=exp(-5*switchind(y2(1:)))*h/f1<cvglobal
	end function

	function getapriorirej(y1,y2,fa1,fa2,tst) result(val)
		type(tsttype)	:: tst
		real	:: y1(0:k),y2(0:k),fa1,fa2
		logical	:: val
		real	:: Zsum,h,var,var1,var2,cv,cvw
		integer	:: j
		Zsum=y1(0)-y2(0)
		Zsum=Zsum+sum(y1(1:k))-sum(y2(1:k))
		var1=sum(y1(1:k)**2)
		var2=sum(y2(1:k)**2)
		var=1.0+var1+var2
		cvw=1.0/var
		cv=cvw*tst%cv_zt(1)+(1.0-cvw)*tst%cv_zt(2)
		
		Zsum=abs(Zsum/sqrt(var))
		val=Zsum>cv
		if(.not. val) return
		if(tst%stage<=1) return
		
		val=rej_stage1(y1,y2,fa1,tst) .and. rej_stage1(y2,y1,fa2,tst)
	end function

	subroutine setcred
		integer	:: l,lb,i
		integer	:: l1,l2,lb2,i2,ix
		integer(8)	:: code
		integer	:: rcount=0
		type(tsttype)	:: tst

		call setcvs
		call settst(tst)

!$omp parallel do private(lb2,l1,l2,i2,code,ix) reduction(+:rcount)
		do lb=1,nblocks
			lb2=lb+1
			if(lb2>nblocks) lb2=1
			do l1=1,npb
				l2=1
				do i2=1,nbi
					code=0
					do ix=1,nfpc
						code=2*code
						if(getapriorirej(y0s(:,l1,lb),y0s(:,l2,lb2),falt(l1,lb),falt(l2,lb2),tst)) then
							code=code+1
							rcount=rcount+1
						endif
						l2=l2+1
					enddo
					code=code*2**(64-nfpc)
					credflags(i2,l1,lb)=code
				enddo
			enddo
		enddo
		print *,"done with setcred"
		call printtime
		print *,"fraction in cred set", real(nsim*npb-rcount)/(nsim*npb)
		ctest=credflags
	end subroutine
	
end module